using FastGaussQuadrature
ngpx = 3; ngpz = 4
gx0 = gausslegendre(ngpx); gaussdatax = [gx0[1] gx0[2]]
gz0 = gausslegendre(ngpz); gaussdataz = [gz0[1] gz0[2]]
gx = gausslegendre(3)[1];
gz = gausslegendre(3)[1];
# geometrical parameters
wtop = 0.3; wbot = 0.45; wweb = 0.015; ftop = 0.02; fbot = 0.03;
hweb = 1.55; Hc = 0.2; Bc = 2.3
h1 = Hc/2.0; h2 = hweb/2.0 + ftop;
span = 25.0 # beam span
# physical parameters
kucs = 4e8; kwcs = 1e5*kucs
μc = 0.2; E28 = 34000e6
Es = 2.1e11; μs = 0.3;
Ar = 0.0; Er = 2.0e11; zr = 0.0;
q0 = -64.56e3;
time0 = [7.0, 30.0, 50.0, 80.0, 120.0, 170.0, 230.0, 300.0, 365.0];
time0 = logspace(log10(30), log10(25550), 30)

function sublength(i::Int64, n::Int64, X0::Float64, X1::Float64, α::Float64)
   s = X1 - X0;
   if abs(α - 1.0) < 1e-6
      return s/(2*n)
   else
      return (s*(1.0 - α))/(2*(1.0 - α^n)) * α^(i - 1)
   end
end

function subnodes(n::Int64, α::Float64, X0::Float64, X1::Float64)
#   X0~X1共划分2n区间
#   X0~(X0+X1)/2之间，从左到右区间长度为公比 α 的等比数列
#   (X0+X1)~X2之间为其左侧数据的镜像
#   返回结果共有2n+1个x坐标
   xcoord = zeros(2*n + 1)
   xcoord[1] = X0
   for i = 2:(n+1)
      xcoord[i] = xcoord[i - 1] + sublength(i - 1, n, X0, X1, α);
   end
   for i = 1:n
      xcoord[n + 1 + i] = xcoord[n + i] + sublength(n + 1 - i, n, X0, X1, α);
   end
   return xcoord
end
